{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook describes the MATLAB package that implements the leave-out correction of <a href=\"https://eml.berkeley.edu/~pkline/papers/KSS2020.pdf\" target=\"_blank\"> Kline, Saggio and Sølvsten (2020)</a> -- KSS henceforth --  for two-way fixed effects models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Introduction\" data-toc-modified-id=\"Introduction-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Introduction</a></span><ul class=\"toc-item\"><li><span><a href=\"#The-KSS-Correction\" data-toc-modified-id=\"The-KSS-Correction-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>The KSS Correction</a></span><ul class=\"toc-item\"><li><span><a href=\"#The-Plug-in-Estimator\" data-toc-modified-id=\"The-Plug-in-Estimator-1.1.1\"><span class=\"toc-item-num\">1.1.1&nbsp;&nbsp;</span>The Plug-in Estimator</a></span></li><li><span><a href=\"#The-Bias-in-the-Plug-in-Estimator\" data-toc-modified-id=\"The-Bias-in-the-Plug-in-Estimator-1.1.2\"><span class=\"toc-item-num\">1.1.2&nbsp;&nbsp;</span>The Bias in the Plug-in Estimator</a></span></li><li><span><a href=\"#The-Problem-with-&quot;Standard&quot;-Standard-Errors-in-High-Dimensional-Models\" data-toc-modified-id=\"The-Problem-with-&quot;Standard&quot;-Standard-Errors-in-High-Dimensional-Models-1.1.3\"><span class=\"toc-item-num\">1.1.3&nbsp;&nbsp;</span>The Problem with \"Standard\" Standard Errors in High Dimensional Models</a></span></li><li><span><a href=\"#The-Leave-Out-Correction\" data-toc-modified-id=\"The-Leave-Out-Correction-1.1.4\"><span class=\"toc-item-num\">1.1.4&nbsp;&nbsp;</span>The Leave Out Correction</a></span></li></ul></li></ul></li><li><span><a href=\"#Computing-the-KSS-Correction\" data-toc-modified-id=\"Computing-the-KSS-Correction-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Computing the KSS Correction</a></span><ul class=\"toc-item\"><li><span><a href=\"#Setup\" data-toc-modified-id=\"Setup-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Setup</a></span></li><li><span><a href=\"#Importing-the-Data\" data-toc-modified-id=\"Importing-the-Data-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Importing the Data</a></span></li><li><span><a href=\"#Calling-the-Main-Function\" data-toc-modified-id=\"Calling-the-Main-Function-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Calling the Main Function</a></span></li></ul></li><li><span><a href=\"#Interpreting-the-Output\" data-toc-modified-id=\"Interpreting-the-Output-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Interpreting the Output</a></span></li><li><span><a href=\"#What-Does-the-Code-Save?\" data-toc-modified-id=\"What-Does-the-Code-Save?-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>What Does the Code Save?</a></span></li><li><span><a href=\"#Scaling-to-Large-Datasets\" data-toc-modified-id=\"Scaling-to-Large-Datasets-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Scaling to Large Datasets</a></span><ul class=\"toc-item\"><li><span><a href=\"#Computational-Bottleneck\" data-toc-modified-id=\"Computational-Bottleneck-5.1\"><span class=\"toc-item-num\">5.1&nbsp;&nbsp;</span>Computational Bottleneck</a></span></li><li><span><a href=\"#The-JLA-Algorithm\" data-toc-modified-id=\"The-JLA-Algorithm-5.2\"><span class=\"toc-item-num\">5.2&nbsp;&nbsp;</span>The JLA Algorithm</a></span></li></ul></li><li><span><a href=\"#Adding-Controls\" data-toc-modified-id=\"Adding-Controls-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Adding Controls</a></span></li><li><span><a href=\"#Leaving-out-a-Person-Year-Observation-vs.-Leaving-Out-a-Match\" data-toc-modified-id=\"Leaving-out-a-Person-Year-Observation-vs.-Leaving-Out-a-Match-7\"><span class=\"toc-item-num\">7&nbsp;&nbsp;</span>Leaving out a Person-Year Observation vs. Leaving Out a Match</a></span></li><li><span><a href=\"#Variance-of-Person-Effects-when-Leaving-Out-a-Match\" data-toc-modified-id=\"Variance-of-Person-Effects-when-Leaving-Out-a-Match-8\"><span class=\"toc-item-num\">8&nbsp;&nbsp;</span>Variance of Person Effects when Leaving Out a Match</a></span></li><li><span><a href=\"#Regressing-Firm-Fixed-Effects-on-Observables\" data-toc-modified-id=\"Regressing-Firm-Fixed-Effects-on-Observables-9\"><span class=\"toc-item-num\">9&nbsp;&nbsp;</span>Regressing Firm Fixed Effects on Observables</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "Economists often study settings where units possess two or more group memberships, some of which can change over time. A prominent example comes from Abowd, Kramarz, and Margolis (1999) (henceforth AKM) who proposed a panel model of log wage determination that is additive in worker and firm fixed effects. \n",
    "\n",
    "This so-called “two-way” fixed effects or \"AKM\" model takes the form:\n",
    "\n",
    "\\begin{equation}\n",
    "\ty_{gt} =  \\alpha_{{g}} + \\psi_{j({g},t)} + w_{gt}'\\delta +  \\varepsilon_{gt}  \\qquad({g}=1,\\dots,N, \\ t=1,\\dots,T_{g} \\ge 2)\n",
    "\\end{equation}\n",
    "\n",
    "where the function $j(\\cdot ,\\cdot ):\\lbrace 1,\\dots ,N\\rbrace \\times \\lbrace 1,\\dots ,\\max_i T_g \\rbrace \\to \\lbrace 1,\\dots ,J\\rbrace$  assigns a worker $g$ and year $t$ observation to one of $J$ firms. Here $\\alpha_g$ is a person effect, $\\psi_{j(g,t)}$ is a firm effect, $w_{gt}$ is a time-varying covariate, and $\\varepsilon_{gt}$ is a time-varying error. \n",
    "\n",
    "Note that can we simply rewrite the original AKM model as:\n",
    "\n",
    "\\begin{equation}\n",
    "y_i =x_i^{\\prime } \\beta +\\varepsilon_i \\qquad i=1,...,n\n",
    "\\end{equation}\n",
    "\n",
    "where $i$ indexes a particular person-year observation $(g,t)$, $x_i$ is a vector that collects all the worker, firm dummies as well as the time-varying covariates $w_{gt}$ so that $\\beta =(\\alpha ,\\psi ,\\delta )'$ is a $k\\times 1$ vector that collects all the worker, firm fixed effects along with $\\delta$.\n",
    "\n",
    "Interest in AKM models often centers on understanding how much of the variability in log wages is attributable to firms. It is common to summarize the firm contribution to wage inequality using the following two parameters:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{\\psi }^2 =\\frac{1}{n}\\sum_{g=1}^N \\sum_{t=1}^{T_g } {\\left(\\psi_{j\\left(g,t\\right)} -\\bar{\\psi} \\right)}^2 \\qquad \\text{and }\\sigma_{\\alpha ,\\psi } =\\frac{1}{n}\\sum_{g=1}^N \\sum_{t=1}^{T_g } \\left(\\psi_{j\\left(g,t\\right)} -\\bar{\\psi} \\right)\\alpha_g\t\t\t\n",
    "\\end{equation}\n",
    "\n",
    "where $\\bar{\\psi} =\\frac{1}{n}\\sum_{g=1}^N \\sum_{t=1}^{T_g } \\psi_{j(g,t)}$. The variance component $\\sigma_{\\psi }^2$ measures the contribution of firm wage variability to inequality, while the covariance component $\\sigma_{\\alpha ,\\psi }$ measures the additional contribution of systematic sorting of high wage workers to high wage firms. \n",
    "\n",
    "The function `leave_out_KSS` provides unbiased estimates of $\\sigma_{\\psi }^2$ and $\\sigma_{\\alpha ,\\psi }$ as well as an estimate of $\\sigma_{\\alpha }^2 =\\frac{1}{n}\\sum_{g=1}^N \\sum_{t=1}^{T_g } {\\left(\\alpha_g -\\bar{\\alpha} \\right)}^2$ using the leave-out bias correction approach proposed by KSS. \n",
    "\n",
    "## The KSS Correction\n",
    "\n",
    "We now provide some general intuition about the KSS leave-out methodology. A more formal discussion can be found in the ECTA article.\n",
    "\n",
    "### The Plug-in Estimator\n",
    "\n",
    "Suppose the researcher is interested in the variance of firm effects, $\\sigma_{\\psi }^2$. To simplify the exposition, we normalize the firm effects so that their firm-size weighted mean is equal to zero, i.e. $\\bar{\\psi}=0$ and re-write $\\sigma_{\\psi }^2$ as \n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{\\psi }^2 =\\frac{1}{J}\\sum_{j=1}^J s_{j}\\psi^{2}_{j}\n",
    "\\end{equation}\n",
    "\n",
    "where $s_{j}$ gives the total number of worker-year observations associated with firm $j$, i.e. $s_{j}=\\sum_{g=1}^{N}\\sum_{t=1}^{T_{g}}\\mathbf{1}\\{j(g,t)=j\\}$.\n",
    "\n",
    "\n",
    "It is customary to report \"plug-in\" estimates of a given variance component using the corresponding OLS estimate. For instance, the plug-in estimate of the variance of firm effects $\\sigma_{\\psi }^2$ is given by\n",
    "\n",
    "\\begin{equation}\n",
    "\\tilde{\\sigma}_{\\psi}^2=\\frac{1}{J}\\sum_{j=1}^J s_{j}\\hat{\\psi}^{2}_{j}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\hat{\\psi}_{j}$ is the OLS estimate obtained after estimating equation (1) via OLS. \n",
    "\n",
    "### The Bias in the Plug-in Estimator\n",
    "\n",
    "\n",
    "The estimated firm effect, $\\hat{\\psi}_{j}$, represents a noisy estimate of the true firm effect, $\\psi_{j}$. The presence of noise in $\\hat{\\psi}_{j}$ is not an issue when one is interested in $\\psi_{j}$ as the OLS estimator $\\hat{\\psi}_{j}$ is assumed to be unbiased in this context, i.e. $E[\\hat{\\psi}_{j}]=\\psi_{j}$.\n",
    "\n",
    "However, estimation error in $\\hat{\\psi}_{j}$ is going to lead to biases when estimating $\\psi_{j}^{2}$ using its \"plug-in\" analogue $\\hat{\\psi}_{j}^{2}$ since:  \n",
    "\n",
    "\\begin{equation}\n",
    "E[\\hat{\\psi}_{j}^{2}]=E[(\\hat{\\psi}_{j}-\\psi_{j}+\\psi_{j})^2]=\\psi^{2}+\\underbrace{\\mathbb{V}[\\hat{\\psi}_{j}]}_{\\text{Bias}}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathbb{V}[\\hat{\\psi}_{j}]$ is the squared standard error of $\\hat{\\psi}_{j}$. Intuitively, when we take the square of $\\hat{\\psi}_{j}$ we are not only squaring its signal, $\\psi_{j}$, but also the estimation error in each $\\hat{\\psi}_{j}$. The latter is going to introduce a bias when estimating $\\psi^{2}_{j}$.   \n",
    "\n",
    "The same logic applies when analyzing the bias of the plug-in estimator of the variance of firm effects since\n",
    "\n",
    "\\begin{equation}\n",
    "E[\\tilde{\\sigma}_{\\psi}^2]=\\sigma^{2}_{\\psi}+\\underbrace{\\dfrac{1}{J}\\sum_{j=1}^{J}s_{j}{\\mathbb{V}[\\hat{\\psi}_{j}]}}_{\\text{Bias}}\n",
    "\\end{equation}\n",
    "\n",
    "### The Problem with \"Standard\" Standard Errors in High Dimensional Models\n",
    "\n",
    "The formula above shows that the bias of the plug-in estimator of the variance of firm effects is\n",
    "\n",
    "\\begin{equation}\n",
    "E[\\tilde{\\sigma}_{\\psi}^2]-\\sigma^{2}_{\\psi}=\\dfrac{1}{J}\\sum_{j=1}^{J}s_{j}{\\mathbb{V}[\\hat{\\psi}_{j}]}.\n",
    "\\end{equation}\n",
    "\n",
    "Therefore, all that is required for a bias correction is an estimate of the (squared) standard error of each firm effect, ${\\mathbb{V}[\\hat{\\psi}_{j}]}$. Similarly, if our interest is on the variance of person effects, then we would need the standard error on each of the person effects, ${\\mathbb{V}[\\hat{\\alpha}_{i}]}$. If we are interested in the covariance of worker and firm effects, then we would need the covariances in sampling errors between each $\\hat{\\alpha}_{i}$ and $\\hat{\\psi}_{j(i,t)}$. \n",
    "\n",
    "The discussion above highlights that an estimate of the sampling variability of the OLS coefficient vector $\\hat{\\beta}=(\\hat{\\alpha},\\hat{\\psi},\\hat{\\delta})$ is required in order to derive an unbiased estimate of the variance components of the AKM model displayed in equation (2). \n",
    "\n",
    "Recall that the sampling variability in $\\hat{\\beta}$, assuming independence across observations, is given by\n",
    "\n",
    "\\begin{equation}\n",
    "{\\mathbb{{V}}}[\\hat{\\beta}]=\\left(\\sum_{i=1}^{n}x_{i}x_{i}'\\right)^{-1}\\sum_{i=1}^{n}{\\sigma}^{2}_{i}x_{i}x_{i}'\\left(\\sum_{i=1}^{n}x_{i}x_{i}'\\right)^{-1}.\n",
    "\\end{equation}\n",
    "\n",
    "where $\\sigma^{2}_{i}=Var(\\varepsilon_{i})$.\n",
    "\n",
    "One can be tempted to provide an estimate of $\\mathbb{V}[\\hat{\\beta}]$ using heteroskedatic *consistent* (\"HC\" or robust), standard errors.  HC standard-errors are calculated using a simple plug-in estimate of $\\sigma^{2}_{i}$ based on\n",
    "\n",
    "\\begin{equation}\n",
    "\\tilde{\\sigma}^{2}_{i}=(y_{i}-x_{i}'\\hat{\\beta})^2.\n",
    "\\end{equation}\n",
    "\n",
    "so that the HC based estimate of $\\mathbb{V}[\\hat{\\beta}]$ is given by \n",
    "\n",
    "\\begin{equation}\n",
    "\\tilde{\\mathbb{{V}}}[\\hat{\\beta}]=\\left(\\sum_{i=1}^{n}x_{i}x_{i}'\\right)^{-1}\\sum_{i=1}^{n}\\tilde{{\\sigma}}^{2}_{i}x_{i}x_{i}'\\left(\\sum_{i=1}^{n}x_{i}x_{i}'\\right)^{-1}.\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "However, HC standard errors based on $\\tilde{\\sigma}^{2}_{i}$  are inconsistent in any high dimensional model where the number of parameters is proportional to the sample size. This corresponds to our context as we often have a firm effect that is estimated using less than 5 observations or \"moves\" of workers across employers.\n",
    "\n",
    "### The Leave Out Correction\n",
    "\n",
    "KSS provide an heteroskedastic-*unbiased* (HU) estimate of the standard error of any coefficient obtained from a high-dimensional linear model.\n",
    "\n",
    "The KSS HU standard error estimate is based on a leave-out estimate of $\\sigma^{2}_{i}$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{\\sigma}^{2}_{i}=y_{i}(y_{i}-x_{i}'\\hat{\\beta}_{-i}).\n",
    "\\end{equation}\n",
    "\n",
    "where $\\hat{\\beta}_{-i}$ is the OLS estimate of $\\beta$ from equation (2) after leaving observation $i$ out. \n",
    "\n",
    "KSS then replaces $\\sigma^{2}_{i}$ in $\\mathbb{V}[\\hat{\\beta}]$ with its unbiased estimate $\\hat{\\sigma}^{2}_{i}$ to derive a HU estimate of $\\mathbb{V}[\\hat{\\beta}]:$\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{\\mathbb{{V}}}[\\hat{\\beta}]=\\left(\\sum_{i=1}^{n}x_{i}x_{i}'\\right)^{-1}\\sum_{i=1}^{n}\\hat{{\\sigma}}^{2}_{i}x_{i}x_{i}'\\left(\\sum_{i=1}^{n}x_{i}x_{i}'\\right)^{-1}.\n",
    "\\end{equation}\n",
    "\n",
    "Going back to the example of the variance of the firm effects, we can extract from $\\hat{\\mathbb{V}}[\\hat{\\beta}]$ the corresponding squared standard error of each firm effect, $\\mathbb{V}^{HU}[\\hat{\\psi}_{j}]$. We can then use the latter to bias correct the corresponding estimate of the variance of the firm effects as follows:\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{\\sigma}^{2}_{\\psi}=\\tilde{\\sigma}^{2}_{\\psi}-\\dfrac{1}{J}\\sum_{j=1}^{J}s_{j}{\\mathbb{V}[\\hat{\\psi}_{j}]}.\n",
    "\\end{equation}\n",
    "\n",
    "The MATLAB function `leave_out_KSS` is going to print the bias corrected variance of firm effects, $\\hat{\\sigma}^{2}_{\\psi}$, as well as the bias-corrected covariance of worker and firm effects and variance of person effects. \n",
    "\n",
    "# Computing the KSS Correction\n",
    "We now demonstrate how one can implement the KSS correction in two-way models using MATLAB and the function `leave_out_KSS`. We continue to work with a simple example based on an AKM model. In what follows, we use the words \"workers\" and \"firms\" when describing the procedure for simplicity but the the function `leave_out_KSS` can be applied to any two-way models (e.g. patients and doctors, students and teachers, strata and treatment arms).\n",
    "\n",
    "## Setup\n",
    "We begin with some auxiliary lines of code that define the relevant paths, call the <a href=\"http://www.cs.cmu.edu/~jkoutis/cmg.html\" target=\"_blank\"> CMG</a> package developed by Yiannis Koutis and set-up the parallel environment within MATLAB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting parallel pool (parpool) using the 'local' profile ...\n",
      "Connected to the parallel pool (number of workers: 6).\n"
     ]
    }
   ],
   "source": [
    "%Setup Paths and Install CMG\n",
    "clc\n",
    "clear\n",
    "cd '/Users/raffaelesaggio/Dropbox/LeaveOutTwoWay'\n",
    "path(path,'codes'); %this contains the main LeaveOut Routines.\n",
    "path(path,'CMG'); % CMG package http://www.cs.cmu.edu/~jkoutis/cmg.html\n",
    "[result,output] = evalc('installCMG(1)'); %installs CMG routine (silently)\n",
    "delete(gcp(\"nocreate\")) %clear parallel envir.\n",
    "c = parcluster('local');  %tell me # of available cores\n",
    "nw = c.NumWorkers; %tell me # of available cores\n",
    "pool=parpool(nw,'IdleTimeout', Inf); %all cores will be assigned to Matlab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing the Data\n",
    "The Github Repo contains a matched employer-employee testing data where we observe the identity of the worker, the identity of the firm employing a given worker, the year in which the match is observed (either 1999 or 2001) and the associated log wage.\n",
    "\n",
    "*Important!:* the original data must be sorted by individual identifiers (id). For instance, one can see that the testing data is sorted by individual identifiers (and year, using `xtset id year` in Stata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Import Data\n",
    "namesrc='data/test.csv'; %path to original testing data\n",
    "data=importdata(namesrc); %import data\n",
    "id=data(:,1); %worker identifiers\n",
    "firmid=data(:,2); %firm identifiers\n",
    "y=data(:,4); % outcome variable\n",
    "clear data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calling the Main Function\n",
    "\n",
    "The function `leave_out_KSS` relies on three mandatory inputs: `(y,id,firmid)`. We can obtain an unbiased variance decomposition of the associated AKM model by simply calling "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Running KSS Correction with the following options\n",
      "Leave Out Strategy: Leave match out\n",
      "Algorithm for Computation of Statistical Leverages: JLA with 200 simulations.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 1\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Info on the leave one out connected set:\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "mean wage: 4.7636\n",
      "variance of wage: 0.1245\n",
      "# of Movers: 6414\n",
      "# of Firms: 1684\n",
      "# of Person Year Observations: 56044\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 2\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Calculating the statistical leverages of the AKM model...\n",
      "Running JLA Algorithm...\n",
      "Done!\n",
      "Elapsed time is 5.602229 seconds.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 3\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "PLUG-IN ESTIMATES (BIASED)\n",
      "Variance of Firm Effects: 0.019821\n",
      "Covariance of Firm, Person Effects: -0.0039091\n",
      "Variance of Person Effects: 0.10354\n",
      "Correlation of Firm, Person Effects: -0.08629\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "BIAS CORRECTED ESTIMATES\n",
      "Variance of Firm Effects: 0.010218\n",
      "Covariance of Firm and Person Effects: 0.0047795\n",
      "Variance of Person Effects: 0.085005\n",
      "Correlation of Firm, Person Effects: 0.16217\n"
     ]
    }
   ],
   "source": [
    "%% Run KSS!\n",
    "[sigma2_psi,sigma_psi_alpha,sigma2_alpha]  = leave_out_KSS(y,id,firmid);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpreting the Output\n",
    "\n",
    "The code `leave_out_KSS` and associated output is composed by three sections. \n",
    "\n",
    "**Section 1**: Here we provide info on leave-out connected set. This is the largest connected set of firms that remains connected after removing any worker from the associated graph, see Lemma 1 and the  Computational Appendix of Kline, Saggio and Sølvsten (2020) for details. The code provides some summary statistics (e.g. # of movers, # of firms, mean and variance of the outcome, etc) for the  leave-out connected set.\n",
    "\n",
    "**Section 2**: After printing the summary statistics, the code computes the statistical leverages of the design, denoted as $P_{ii}$. Computation of $\\{P_{ii}\\}_{i=1}^{n}$  represents the main computational bottleneck of the routine. \n",
    "\n",
    "**Section 3**: The code then enters its third, and final, stage where the main results are printed. The code starts by reporting the --- biased --- estimates of the variance components that result from the \"plug-in\" approach of treating OLS estimates as measured without error. Finally, the code prints the bias corrected variance of firm effects and the covariance of worker and firm effects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What Does the Code Save?\n",
    "\n",
    "`leave_out_KSS` saves three scalars: the variance of firm effects (`sigma2_psi` in [4]), the covariance of worker and firm effects (`sigma_psi_alpha`) and the variance of person effects (`sigma2_alpha`). \n",
    "\n",
    "`leave_out_KSS` also saves on disk one .csv file. This .csv contains information on the leave-out connected set. This file has 4 columns. First column reports the outcome variable, second and third columns the worker and the firm identifiers (as originally inputted by the user). The fourth column reports the statistical leverages of the AKM model. If the code is reporting a leave-out correction at the match-level, the .csv will be collapsed at the match level. By default, the .csv file is going to be saved in the main directory under the name `leave_out_estimates`. The user can specify an alternative path using the option `filename` when calling `leave_out_KSS`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scaling to Large Datasets\n",
    "\n",
    "`leave_out_KSS` can be used on extremely large datasets. The code uses a variant of the random projection method, denoted as the Johnson–\n",
    "Lindenstrauss Approximation (JLA henceforth) algorithm in KSS for its connection to the work of <a href=\"https://www.researchgate.net/profile/William_Johnson16/publication/235008656_Extensions_of_Lipschitz_maps_into_a_Hilbert_space/links/55e9abf908aeb65162649527/Extensions-of-Lipschitz-maps-into-a-Hilbert-space.pdf\" target=\"_blank\"> Johnson and Lindenstrauss (1984)</a>, see also <a href=\"https://www.sciencedirect.com/science/article/pii/S0022000003000254\" target=\"_blank\"> Achlioptas (2003)</a>. We now discuss briefly the main computational bottleneck of the procedure and the JLA algorithm. \n",
    "\n",
    "## Computational Bottleneck\n",
    "\n",
    "Recall from the discussion in Section 1 that the KSS leave-out bias correction procedure relies on leave-out estimates of $\\sigma^{2}_{i}$\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{\\sigma}^{2}_{i}=y_{i}(y_{i}-x_{i}'\\hat{\\beta}_{-i})\n",
    "\\end{equation}\n",
    "\n",
    "where $\\hat{\\beta}_{-i}$ is the OLS estimate of $\\beta$ from the AKM model in equation (2) after leaving observation $i$ out.\n",
    "\n",
    "Clearly, re-estimating $\\hat{\\beta}_{-i}$ by leaving a particular observation $i$ for $n$ times, is infeasible computationally. Luckily, one can re-write $\\hat{\\sigma}^{2}_{i}$ as\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{\\sigma}^{2}_{i}=y_{i}\\frac{(y_{i}-x_{i}'\\hat{\\beta})}{1-P_{ii}}\n",
    "\\end{equation}\n",
    "\n",
    "where $P_{ii}$ measures the influence or leverage of observation $i$, i.e. $P_{ii} =x_i^{\\prime } S_{xx}^{-1} x_i$. The expression above highlights that all that is needed for computation of $\\hat{\\sigma}^{2}_{i}$ are the $n$ statistical leverages $\\{P_{ii}\\}_{i=1}^{n}$. \n",
    "However, exact computation of $P_{ii}$ might still remain prohibitive when $n$ is in order of tens of millions or higher. \n",
    "\n",
    "## The JLA Algorithm\n",
    "The JLA algorithm introduced by KSS provides a stochastic approximation to $\\{P_{ii}\\}_{i=1}^{n}$ using the random projection ideas developed by <a href=\"https://www.researchgate.net/profile/William_Johnson16/publication/235008656_Extensions_of_Lipschitz_maps_into_a_Hilbert_space/links/55e9abf908aeb65162649527/Extensions-of-Lipschitz-maps-into-a-Hilbert-space.pdf\" target=\"_blank\"> Johnson and Lindenstrauss (1984)</a>. We defer the reader to the <a href=\"https://www.dropbox.com/s/ycvls8pbtxewj06/DataComputationAppendix.pdf?dl=1\" target=\"_blank\"> Computational Appendix of KSS</a>  for details.\n",
    "\n",
    "The number of simulations underlying the JLA algorithm is governed by the input `simulations_JLA` (which is defined as $p$ in the computational appendix of KSS). Intuitively, more simulations imply a higher accuracy -- but higher computation time --- when estimating $\\lbrace P_{ii} ,B_{ii} \\rbrace_{i=1}^n$. \n",
    "\n",
    "**Note:** The user might want to pre-specify a random number generator seed to ensure replicability when calling the function  `leave_out_KSS`. \n",
    "\n",
    "We now demonstrate the performance of the code on a large dataset\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Running KSS Correction with the following options\n",
      "Leave Out Strategy: Leave match out\n",
      "Algorithm for Computation of Statistical Leverages: JLA with 50 simulations.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 1\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Info on the leave one out connected set:\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "mean wage: 4.7304\n",
      "variance of wage: 0.16248\n",
      "# of Movers: 916632\n",
      "# of Firms: 165360\n",
      "# of Person Year Observations: 13860616\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 2\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Calculating the statistical leverages of the AKM model...\n",
      "Running JLA Algorithm...\n",
      "Done!\n",
      "Elapsed time is 251.168051 seconds.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 3\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "PLUG-IN ESTIMATES (BIASED)\n",
      "Variance of Firm Effects: 0.039448\n",
      "Covariance of Firm, Person Effects: 0.0084313\n",
      "Variance of Person Effects: 0.080329\n",
      "Correlation of Firm, Person Effects: 0.14978\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "BIAS CORRECTED ESTIMATES\n",
      "Variance of Firm Effects: 0.030371\n",
      "Covariance of Firm and Person Effects: 0.014639\n",
      "Variance of Person Effects: 0.048803\n",
      "Correlation of Firm, Person Effects: 0.38024\n"
     ]
    }
   ],
   "source": [
    "%% Running KSS on a large dataset          \n",
    "websave('large_fake.csv', 'https://www.dropbox.com/s/ny5tef29ij7ran2/large_fake_data.csv?dl=1'); %downloads and saves to disk a fake, large matched employer employee data\n",
    "namesrc='large_fake.csv'; %path to the large data\n",
    "data=importdata(namesrc); %import data \n",
    "id=data(:,1); %worker identifiers\n",
    "firmid=data(:,2); %firm identifiers\n",
    "y=data(:,4); % outcome variable\n",
    "clear data\n",
    "delete('large_fake.csv'); %delete original .csv data from disk\n",
    "\n",
    "%Run Leave Out Correction (50 simulations) \n",
    "type_of_algorithm='JLA'; %run random projection algorithm\n",
    "simulations_JLA=50;\n",
    "[sigma2_psi,sigma_psi_alpha,sigma2_alpha]  = leave_out_KSS(y,id,firmid,[],[],type_of_algorithm,simulations_JLA);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " We can see from the output that the leave-out connected set has almost 14 million person-year observations. The code is able to complete in around 4 minutes (on a 2020 Macbook pro with 6 cores and 16GB of RAM). \n",
    "\n",
    "The computational appendix in KSS shows that the JLA algorithm can cut computation time by a factor of 100 while introducing an approximation error of roughly $10^{-4}$.\n",
    "\n",
    "The current code uses an improved estimator of both $P_{ii}$ and $M_{ii}=1-P_{ii}$  which are both guaranteed to lie in $[0;1]$. These improved estimators are then combined to derive an unbiased JLA estimator of a given variance component provided that $\\frac{n}{p^{4}}=o(1)$, see <a href=\"https://www.dropbox.com/s/i28yvzae2tnp2tl/improved_JLA.pdf?dl=1\" target=\"_blank\">this document</a> for details.  \n",
    "\n",
    "One can check the stability of the estimates for different values of `simulations_JLA`. For instance, if we double `simulations_JLA` from 50 to 100 and run the code again on the same data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Running KSS Correction with the following options\n",
      "Leave Out Strategy: Leave match out\n",
      "Algorithm for Computation of Statistical Leverages: JLA with 100 simulations.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 1\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Info on the leave one out connected set:\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "mean wage: 4.7304\n",
      "variance of wage: 0.16248\n",
      "# of Movers: 916632\n",
      "# of Firms: 165360\n",
      "# of Person Year Observations: 13860616\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 2\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Calculating the statistical leverages of the AKM model...\n",
      "Running JLA Algorithm...\n",
      "Done!\n",
      "Elapsed time is 458.093692 seconds.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 3\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "PLUG-IN ESTIMATES (BIASED)\n",
      "Variance of Firm Effects: 0.039448\n",
      "Covariance of Firm, Person Effects: 0.0084313\n",
      "Variance of Person Effects: 0.080329\n",
      "Correlation of Firm, Person Effects: 0.14978\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "BIAS CORRECTED ESTIMATES\n",
      "Variance of Firm Effects: 0.030382\n",
      "Covariance of Firm and Person Effects: 0.014598\n",
      "Variance of Person Effects: 0.048738\n",
      "Correlation of Firm, Person Effects: 0.37937\n"
     ]
    }
   ],
   "source": [
    "%% Compute estimates while doubling number of simulations\n",
    "simulations_JLA=100;\n",
    "[sigma2_psi,sigma_psi_alpha,sigma2_alpha]  = leave_out_KSS(y,id,firmid,[],[],type_of_algorithm,simulations_JLA); %check stability of variance components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We obtain virtually the same variance components as when `simulations_JLA`=50 while significantly increasing the computational time! If the user does not specify a value for `simulations_JLA`,the code defaults to `simulations_JLA`=200.\n",
    "\n",
    "We conclude this section by noting that the user can also calculate an exact version of $\\lbrace P_{ii} \\rbrace_{i=1}^n$. This can be done by setting the option `type_of_algorithm` to `exact`. \n",
    "\n",
    "**Warning:** Calling the option `exact` in large datasets can be very time consuming! We now load again the original, smaller, testing data and then compare the exact and JLA based estimates of the variance components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Running KSS Correction with the following options\n",
      "Leave Out Strategy: Leave match out\n",
      "Algorithm for Computation of Statistical Leverages: Exact\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 1\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Info on the leave one out connected set:\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "mean wage: 4.7636\n",
      "variance of wage: 0.1245\n",
      "# of Movers: 6414\n",
      "# of Firms: 1684\n",
      "# of Person Year Observations: 56044\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 2\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Calculating the statistical leverages of the AKM model...\n",
      "Running Exact Algorithm...\n",
      "Done!\n",
      "Elapsed time is 162.242224 seconds.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 3\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "PLUG-IN ESTIMATES (BIASED)\n",
      "Variance of Firm Effects: 0.019821\n",
      "Covariance of Firm, Person Effects: -0.0039091\n",
      "Variance of Person Effects: 0.10354\n",
      "Correlation of Firm, Person Effects: -0.08629\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "BIAS CORRECTED ESTIMATES\n",
      "Variance of Firm Effects: 0.010289\n",
      "Covariance of Firm and Person Effects: 0.0046293\n",
      "Variance of Person Effects: 0.085204\n",
      "Correlation of Firm, Person Effects: 0.15635\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Running KSS Correction with the following options\n",
      "Leave Out Strategy: Leave match out\n",
      "Algorithm for Computation of Statistical Leverages: JLA with 100 simulations.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 1\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Info on the leave one out connected set:\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "mean wage: 4.7636\n",
      "variance of wage: 0.1245\n",
      "# of Movers: 6414\n",
      "# of Firms: 1684\n",
      "# of Person Year Observations: 56044\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 2\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Calculating the statistical leverages of the AKM model...\n",
      "Running JLA Algorithm...\n",
      "Done!\n",
      "Elapsed time is 2.434287 seconds.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 3\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "PLUG-IN ESTIMATES (BIASED)\n",
      "Variance of Firm Effects: 0.019821\n",
      "Covariance of Firm, Person Effects: -0.0039091\n",
      "Variance of Person Effects: 0.10354\n",
      "Correlation of Firm, Person Effects: -0.08629\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "BIAS CORRECTED ESTIMATES\n",
      "Variance of Firm Effects: 0.01044\n",
      "Covariance of Firm and Person Effects: 0.0044957\n",
      "Variance of Person Effects: 0.085326\n",
      "Correlation of Firm, Person Effects: 0.15063\n"
     ]
    }
   ],
   "source": [
    "%% Compare Exact vs. JLA Estimates\n",
    "namesrc='data/test.csv'; %path to original testing data\n",
    "data=importdata(namesrc); %import data\n",
    "id=data(:,1); %worker identifiers\n",
    "firmid=data(:,2); %firm identifiers\n",
    "y=data(:,4); % outcome variable\n",
    "clear data\n",
    "\n",
    "%Run Leave Out Correction with exact\n",
    "type_of_algorithm='exact'; %run random projection algorithm;\n",
    "[sigma2_psi,sigma_psi_alpha,sigma2_alpha]  = leave_out_KSS(y,id,firmid,[],[],type_of_algorithm);\n",
    "\n",
    "%Run Leave Out Correction with JLA\n",
    "simulations_JLA=100;\n",
    "type_of_algorithm='JLA'; %run random projection algorithm;\n",
    "[sigma2_psi,sigma_psi_alpha,sigma2_alpha]  = leave_out_KSS(y,id,firmid,[],[],type_of_algorithm,simulations_JLA);\n",
    "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variance components estimated via JLA are extremely close to the `exact` ones but only take a fraction of the time to compute. If the input data has more than 10,000 obs, the code defaults to using the JLA algorithm unless the user specifies type_of_algorithm to \"exact\". "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding Controls\n",
    "\n",
    "We have demonstrated the functioning of `leave_out_KSS` using a simple AKM model with no controls ($w_{gt}=0$). It is easy to add a matrix of controls to the routine. Suppose for instance that we want to add year fixed effects to the original AKM model. This can be done as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Running KSS Correction with the following options\n",
      "Leave Out Strategy: Leave match out\n",
      "Algorithm for Computation of Statistical Leverages: JLA with 200 simulations.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 1\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Info on the leave one out connected set:\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "mean wage: 4.7636\n",
      "variance of wage: 0.1245\n",
      "# of Movers: 6414\n",
      "# of Firms: 1684\n",
      "# of Person Year Observations: 56044\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 2\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "pcg converged at iteration 58 to a solution with relative residual 8.7e-11.\n",
      "Calculating the statistical leverages of the AKM model...\n",
      "Running JLA Algorithm...\n",
      "Done!\n",
      "Elapsed time is 4.271175 seconds.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 3\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "PLUG-IN ESTIMATES (BIASED)\n",
      "Variance of Firm Effects: 0.019479\n",
      "Covariance of Firm, Person Effects: -0.004008\n",
      "Variance of Person Effects: 0.10404\n",
      "Correlation of Firm, Person Effects: -0.089031\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "BIAS CORRECTED ESTIMATES\n",
      "Variance of Firm Effects: 0.0097812\n",
      "Covariance of Firm and Person Effects: 0.0046239\n",
      "Variance of Person Effects: 0.08578\n",
      "Correlation of Firm, Person Effects: 0.15963\n"
     ]
    }
   ],
   "source": [
    "%% How to add controls\n",
    "namesrc='data/test.csv'; %path to original testing data\n",
    "data=importdata(namesrc); %import data\n",
    "id=data(:,1); %worker identifiers\n",
    "firmid=data(:,2); %firm identifiers\n",
    "year=data(:,3); %year identifier\n",
    "y=data(:,4); % outcome variable\n",
    "clear data\n",
    "\n",
    "%Specify year fixed effects as controls\n",
    "[~,~,controls] = unique(year);\n",
    "controls \t   = sparse((1:size(y,1))',controls',1,size(y,1),max(controls));\n",
    "controls       = controls(:,1:end-1); %to avoid collinearity issues, omit last year fixed effects.\n",
    "\n",
    "%Call KSS with matrix of controls\n",
    "[sigma2_psi,sigma_psi_alpha,sigma2_alpha]  = leave_out_KSS(y,id,firmid,controls);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When controls are specified, the code proceeds by partialling them out. That is, it first estimates the AKM model in the leave-out connected set\n",
    "\n",
    "\\begin{equation}\n",
    "y_{gt}=\\alpha_{g}+\\psi_{j(g,t)}+w_{gt}'\\delta+\\varepsilon_{gt}\n",
    "\\end{equation}\n",
    "\n",
    "from which we obtain $\\hat{\\delta}$. We then work with a residualized model where the outcome variable is now defined as $y_{gt}^{new}=y_{gt}-w_{gt}'\\hat{\\delta}$  and project this residualized outcome on worker and firm indicators and report the associated (bias-corrected) variance components. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Leaving out a Person-Year Observation vs. Leaving Out a Match\n",
    "\n",
    "By default, the code reports leave-out corrections for the variance of firm effects and the covariance of firm and worker effects that are robust to unrestricted heteroskedasticity and serial correlation of the error term  within a given match (unique combination of worker and firm identifier), see Remark 3 of KSS. We discuss the interpretation of the leave-out corrected variance of person effects when leaving a match out here. \n",
    "\n",
    "The user can specify the function to run the KSS correction when leaving only an observation out using the option `leave_out_level`. When leaving a person-year observation out, the resulting KSS variance components are robust to unrestricted heteroskedasticity but not serial correlation within match. Below we demonstrate how to compute KSS adjusted variance components when leaving a single (person-year) observation out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Running KSS Correction with the following options\n",
      "Leave Out Strategy: Leave person-year observation out\n",
      "Algorithm for Computation of Statistical Leverages: JLA with 200 simulations.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 1\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Info on the leave one out connected set:\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "mean wage: 4.7636\n",
      "variance of wage: 0.1245\n",
      "# of Movers: 6414\n",
      "# of Firms: 1684\n",
      "# of Person Year Observations: 56044\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 2\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Calculating the statistical leverages of the AKM model...\n",
      "Running JLA Algorithm...\n",
      "Done!\n",
      "Elapsed time is 4.700531 seconds.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 3\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "PLUG-IN ESTIMATES (BIASED)\n",
      "Variance of Firm Effects: 0.019821\n",
      "Covariance of Firm, Person Effects: -0.0039091\n",
      "Variance of Person Effects: 0.10354\n",
      "Correlation of Firm, Person Effects: -0.08629\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "BIAS CORRECTED ESTIMATES\n",
      "Variance of Firm Effects: 0.010306\n",
      "Covariance of Firm and Person Effects: 0.0046087\n",
      "Variance of Person Effects: 0.085253\n",
      "Correlation of Firm, Person Effects: 0.15548\n"
     ]
    }
   ],
   "source": [
    "%% Leaving out a Person-Year Observation vs. Leaving Out a Match \n",
    "\n",
    "leave_out_level='obs'; %leave a single person-year observation out\n",
    "[sigma2_psi,sigma_psi_alpha,sigma2_alpha]  = leave_out_KSS(y,id,firmid,[],leave_out_level);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When $T=2$ (i.e the underlying matched employer-employee data spans only two years), as in this example, it turns out that the KSS adjusted variance of firm effects and covariance of firm and worker effects is robust to any arbitrary correlation between $\\varepsilon_{g2}$ and $\\varepsilon_{g1}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variance of Person Effects when Leaving Out a Match\n",
    "\n",
    "By leaving a match-out, we can bias correct for the variance of firm effects and covariance of worker and firm effects while allowing for unrestricted hetoreskedasticity and serial correlation of the error term $\\varepsilon_{gt}$ within each worker-firm match.\n",
    "\n",
    "However, the person effects, $\\alpha_{g}$, of \"stayers\" --- workers that never leave a particular firm ---  are not leave-match-out estimable. This implies that we cannot compute an unbiased estimate of $\\Omega_{g}=Var(\\varepsilon_{g1},...,\\varepsilon_{gT_{g}})$ for stayers. An estimate of $\\Omega_{g}$ for both stayers and movers is required in order to provide a bias correction for the variance of person effects, see Section 1 and Remark 3 in KSS.\n",
    "\n",
    "The current implementation of the code provides an estimate of $\\Omega_{g}$ for stayers that is robust to unrestricted heteroskedasticity while using an estimate of $\\Omega_{g}$ for movers that is robust to both unrestricted heteroskedasticity and serial correlation within a match. This allows us to compute an estimate of the variance of person effects that represents an upper bound estimate on the variance of person effects (computed across both stayers and movers).   \n",
    "\n",
    "There are several alternatives that the user can explore: \n",
    "\n",
    "\n",
    "1. Estimate a variance decomposition in a sample of movers only: For movers, it is possible to estimate a leave-out bias corrected variance of person effects that is robust to both unrestricted heteroskedasticity and serial correlation in the error term of the AKM model within a given match. Therefore, one can provide an unbiased variance decomposition of all the three components of the two-way fixed effects model by simply feeding to the function `leave_out_KSS` a movers-only sample.\n",
    "\n",
    "\n",
    "\n",
    "2. Drop adjacent wage observations for stayers: Under the assumption that the errors are serially independent after m periods, it suffices to keep every $mth$ stayer observation and apply the leave person-year out estimator.  For example, if $m=2$ and we have a balanced panel with $T=5$, we can restore independence of the errors in the stayer sample by keeping any of the following pairs of stayer time periods: (1,4), (2,5), (1,5). One can choose from the available pairs randomly for each stayer with equal probability.    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regressing Firm Fixed Effects on Observables\n",
    "\n",
    "It is common in empirical applications to regress the fixed effects estimated from the two-way model on some observables characteristics. Using the AKM model again as our leading example, suppose we are interested in the linear projection of the firm effects, $\\psi_{gt}$, on some observables, $Z_{gt}$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\psi_{j(g,t)}=Z_{gt}'\\gamma+e_{gt}\n",
    "\\end{equation}\n",
    "\n",
    "Typical practice is to estimate $\\gamma$ using a simple regression where the estimated firm effects, $\\hat{\\psi}_{j(g,t)}$, are regressed on $Z_{gt}$\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{\\gamma}=\\left(\\sum_{g,t}Z_{gt}Z_{gt}'.\n",
    "\\right)^{-1}Z_{gt}\\hat{\\psi}_{gt}\n",
    "\\end{equation}\n",
    "\n",
    "KSS show that inference on $\\hat{\\gamma}$ needs to be adjusted. The reason is that the firm fixed effects $\\hat{\\psi}_{j(g,t)}$ are all potentially correlated with one another. \n",
    "\n",
    "To see this, suppose that we have a simple AKM model with only two time periods, set $w_{gt}=0$, and take first differences $\\Delta y_{g}\\equiv y_{g2}-y_{g1}$  to eliminate the worker fixed effects so that the AKM model becomes\n",
    "\\begin{equation}\n",
    "\\Delta y_{g}=\\Delta f_{g}'\\psi+\\varepsilon_{g}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\Delta f_{g}=f_{g,2}-f_{g,1}$ and $f_{gt}=\\{\\mathbf{1}_{j(g,t)=1},..,\\mathbf{1}_{j(g,t)=J}\\}$ is the vector containing the firm dummies. In this model,\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{\\psi}=\\psi+\\underbrace{\\sum_{g=1}^{N}(\\Delta f_{g}\\Delta f_{g}')^{-1}\\Delta f_{g}\\varepsilon_{g}}_{\\text{Correlated Noise}}\n",
    "\\end{equation}\n",
    "\n",
    "Notice how the dependence in the vector of estimated firm fixed effects, $\\hat{\\psi}$, is induced by the design $\\sum_{g=1}^{N}(\\Delta f_{g}\\Delta f_{g}')^{-1}$.  Intuitively, conducting inference on $\\hat{\\gamma}$ is challenging as we have to provide standard errors in a context where all of our underlying observations are potentially correlated with one another. As shown in Table 3 of KSS, ignoring this correlation can easily lead standard errors to be underestimated by an order of magnitude in practice.\n",
    "\n",
    "The package provides the correct standard errors on $\\hat{\\gamma}$ using the function `lincom_KSS` which is designed like the Stata function `lincom` and therefore works as a post-estimation command. We demonstrate the functioning of `lincom_KSS` with an example.\n",
    "\n",
    "In this example, we are interested in testing whether the difference in person-year weighted mean firm effects between region 1 and region 2 is statistically different from zero. This amounts to running a regression where the dependent variable is the vector of estimated firm effects and the set of observables, $Z_{gt}$ , here is represented by a constant and a dummy for whether the firm of worker $g$ in year $t$ belongs to region 2.\n",
    "\n",
    "The resulting coefficient (and standard error) can be computed by calling the function `leave_out_KSS` specifying that we want to run the `lincom` option and using the region dummy as $Z_{gt}$ (the constant is automatically added by the code)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Running KSS Correction with the following options\n",
      "Leave Out Strategy: Leave match out\n",
      "Algorithm for Computation of Statistical Leverages: JLA with 200 simulations.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 1\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Info on the leave one out connected set:\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "mean wage: 4.7047\n",
      "variance of wage: 0.14653\n",
      "# of Movers: 9972\n",
      "# of Firms: 2974\n",
      "# of Person Year Observations: 89666\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 2\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "Calculating the statistical leverages of the AKM model...\n",
      "Running JLA Algorithm...\n",
      "Done!\n",
      "Elapsed time is 9.324213 seconds.\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "SECTION 3\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "PLUG-IN ESTIMATES (BIASED)\n",
      "Variance of Firm Effects: 0.060695\n",
      "Covariance of Firm, Person Effects: -0.012603\n",
      "Variance of Person Effects: 0.10318\n",
      "Correlation of Firm, Person Effects: -0.15926\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n",
      "BIAS CORRECTED ESTIMATES\n",
      "Variance of Firm Effects: 0.044613\n",
      "Covariance of Firm and Person Effects: 0.0025688\n",
      "Variance of Person Effects: 0.079191\n",
      "Correlation of Firm, Person Effects: 0.043218\n",
      "Regressing the firm effects on observables...\n",
      "pcg converged at iteration 115 to a solution with relative residual 8.6e-11.\n",
      "******************************************\n",
      "******************************************\n",
      "RESULTS ON LINCOM\n",
      "******************************************\n",
      "******************************************\n",
      "Coefficient on Region 2 Dummy:  0.25982\n",
      "Robust \"White\" Standard Error: 0.050155\n",
      "KSS Standard error:  0.088374\n",
      "T-stat: 2.94\n",
      "******************************************\n"
     ]
    }
   ],
   "source": [
    "%Regressing firm effects on observables \n",
    "namesrc='data/lincom.csv'; %testing data for the lincom function\n",
    "data=importdata(namesrc);\n",
    "id=data(:,1); \n",
    "firmid=data(:,2);\n",
    "y=data(:,5);\n",
    "region=data(:,4); %Region indicator. Value -1 for region 1, Value 1 for region 2;\n",
    "region_dummy=region;\n",
    "region_dummy(region_dummy==-1)=0; %Make it a proper dummy variable\n",
    "\n",
    "%Run the KSS correction and \"lincom\"\n",
    "labels_lincom={'Region 2 Dummy'}; %give me the label of the columns of Z.\n",
    "lincom_do=1; %tell the function leave_out_KSS that we want to project the firm effects on some Z.\n",
    "Z=region_dummy; %we're going to project the firm effects on a constant +  the region dummy. Constant automatically added by the code\n",
    "\n",
    "%Ready to call KSS with lincom option!\n",
    "[sigma2_psi,sigma_psi_alpha,sigma2_alpha] = leave_out_KSS(y,id,firmid,[],[],[],[],lincom_do,Z,labels_lincom);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the output above (make sure to scroll until the end) that the difference in person-year weighted mean firm effects between the two regions is equal to 0.26. The traditional \"heteroskedatic-consistent\" or  \"robust\" standard errors on this coefficient is around 0.05 while the KSS-adjusted standard error is roughly twice as large."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "matlab",
   "version": "0.16.11"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
